﻿#pragma once

#include <learn/geom/vec.h>
#include <learn/math/math.h>

namespace xi
{

namespace geom
{

// 使用 const 指针是为了能够在 const 成员函数中使用，因为 const 成员函数不能修改成员变量。
// 使用 int 而非 std::size_t 是为了能够避免下标运算时的类型转换问题。

/**
 * @brief 计算幂基曲线的点
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] n 控制点数
 * @param[in] u
 * @return EXT_SYMBOL
 */
template <class T> EXT_SYMBOL T horner1(const T *cpts, int n, double u);

/**
 * @brief 计算幂基曲面上的点
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] n u 方向控制点数
 * @param[in] m v 方向控制点数
 * @param[in] u
 * @param[in] v
 * @return EXT_SYMBOL
 */
template <class T> EXT_SYMBOL T horner2(const T *cpts, int n, int m, double u, double v);

/**
 * @brief 计算伯恩斯坦多项式 B^n_i(u) = u^i(1-u)^(n-i)
 *
 * @param[in] n
 * @param[in] i
 * @param[in] u
 * @return EXT_SYMBOL
 */
EXT_SYMBOL double bernstein(int n, int i, double u);

/**
 * @brief 计算所有伯恩斯坦多项式 B^n_i(u) = u^i(1-u)^(n-i)
 *
 * @param[in] n
 * @param[in] u
 * @return EXT_SYMBOL
 */
EXT_SYMBOL std::vector<double> all_bernstein(int n, double u);

/**
 * @brief 计算贝塞尔曲线上的点
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] n 控制点数
 * @param[in] u 参数
 * @return EXT_SYMBOL
 */
template <class T> EXT_SYMBOL T bezier_point(const T *cpts, int n, double u);

/**
 * @brief 使用 de Casteljau 算法计算贝塞尔曲线上的点
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] n 控制点数
 * @param[in] u 参数
 * @return EXT_SYMBOL
 */
template <class T> EXT_SYMBOL T de_casteljau1(const T *cpts, int n, double u);

/**
 * @brief 使用 de Casteljau 算法计算贝塞尔曲面上的点
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] n u 方向控制点数
 * @param[in] m v 方向控制点数
 * @param[in] u 参数
 * @param[in] v 参数
 * @return EXT_SYMBOL
 */
template <class T> EXT_SYMBOL T de_casteljau2(const T *cpts, int n, int m, double u, double v);

/**
 * @brief Bezier 曲线升次
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] n 控制点数
 * @param[in] t 升次次数
 */
template <class T> EXT_SYMBOL std::vector<T> bezier_elevate(const T *cpts, int n, int t);

/**
 * @brief Bezier 曲线降次
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] n 控制点数
 * @param[out] Q 降次后的控制点
 * @return EXT_SYMBOL
 */
template <class T> EXT_SYMBOL double bezier_reduce(const T *cpts, int n, T *Q);

/**
 * @brief 幂基曲线转 Bezier 曲线
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] n 控制点数
 * @param[out] R Bezier 曲线控制点
 * @return EXT_SYMBOL
 */
template <class T> EXT_SYMBOL void power_to_bezier(const T *cpts, int n, T *R);

/**
 * @brief Bezier 转幂基矩阵
 *
 * @param[in] p Bezier 阶数
 * @param[in] M 幂基转 Bezier 矩阵
 * @return EXT_SYMBOL
 */
EXT_SYMBOL void bezier_to_power_matrix(int p, double *M);

/**
 * @brief 幂基转 Bezier 矩阵
 *
 * @param[in] p 次数
 * @param[in] M 幂基转 Bezier 矩阵
 * @param[out] MI 幂基转 Bezier 矩阵的逆矩阵
 * @return EXT_SYMBOL
 */
EXT_SYMBOL void power_to_bezier_matrix(int p, const double *M, double *MI);

/**
 * @brief 搜索参数 u 在 p 次 knots 序列中的位置
 *
 * @param[in] knots 节点向量
 * @param[in] n 节点向量长度
 * @param[in] p 次数
 * @param[in] u
 * @return EXT_SYMBOL
 */
EXT_SYMBOL int find_span(const double *knots, int n, int p, double u);

/**
 * @brief 搜索参数 u 在 p 次 knots 序列中的重数
 *
 * @param[in] knots 节点向量
 * @param[in] i u 在 knots[i] 和 knots[i+1] 之间
 * @param[in] p 次数
 * @param[in] u
 * @return EXT_SYMBOL
 */
EXT_SYMBOL int find_mult_span(const double *knots, int i, int p, double u);

/**
 * @brief 计算所有非零基函数 N_{i,p}(u) : 下标范围 basis[j] = N_{i-p+j}(u)
 *
 * @param[in] knots 节点向量
 * @param[in] p 次数
 * @param[in] i 基函数编号
 * @param[in] u
 * @return EXT_SYMBOL
 */
EXT_SYMBOL std::vector<double> basis_funcs(const double *knots, int p, int i, double u);

/**
 * @brief 计算所有非零基函数 N_{i,k}(u) k = 0,1,...,p : 下标范围 basis[k][j] = N_{i-k+j,k}(u)
 *
 * @param[in] knots 节点向量
 * @param[in] p 次数
 * @param[in] i 基函数编号
 * @param[in] u
 * @return EXT_SYMBOL
 */
EXT_SYMBOL std::vector<std::vector<double>> all_basis_funcs(const double *knots, int p, int i, double u);

/**
 * @brief 计算所有非零基函数的导数 N_{i,p}^(d)(u) : ders[k][j] 表示 N_{i-p+j}^(k)
 *
 * @param[in] knots 节点向量
 * @param[in] p 次数
 * @param[in] i 基函数编号
 * @param[in] d 导数阶数
 * @param[in] u
 * @return EXT_SYMBOL
 */
EXT_SYMBOL std::vector<std::vector<double>> deriv_basis_funcs(const double *knots, int p, int i, int d, double u);

/**
 * @brief 计算单个基函数 N_{i,p}(u)
 *
 * @param[in] knots 节点向量
 * @param[in] n 节点向量长度
 * @param[in] p 次数
 * @param[in] i 基函数编号
 * @param[in] u
 * @return EXT_SYMBOL
 */
EXT_SYMBOL double one_basis_func(const double *knots, int n, int p, int i, double u);

/**
 * @brief 计算单个基函数不超过 d 阶导数 N_{i,p}^(d)(u) : ders[k] 表示 N_{i,p}^(k)，i 必须是 find_span 找到的下标位置
 *
 * @param[in] knots 节点向量
 * @param[in] n 节点向量长度
 * @param[in] p 次数
 * @param[in] i 基函数编号
 * @param[in] d 导数阶数
 * @param[in] u
 * @return EXT_SYMBOL
 */
EXT_SYMBOL std::vector<double> deriv_one_basis_func(const double *knots, int n, int p, int i, int d, double u);

/**
 * @brief 计算 B 样条曲线上的点
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] knots 节点向量
 * @param[in] n 节点向量长度
 * @param[in] p 次数
 * @param[in] u
 * @return EXT_SYMBOL
 */
template <class T> EXT_SYMBOL T curve_point(const T *cpts, const double *knots, int n, int p, double u);

/**
 * @brief 计算 B 样条曲面上的点
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] uknots u 方向节点向量
 * @param[in] vknots v 方向节点向量
 * @param[in] n u 方向节点向量长度
 * @param[in] m v 方向节点向量长度
 * @param[in] p u 方向次数
 * @param[in] q v 方向次数
 * @param[in] u
 * @param[in] v
 * @return EXT_SYMBOL
 */
template <class T>
EXT_SYMBOL T surface_point(const T *cpts, const double *uknots, const double *vknots, int n, int m, int p, int q,
                           double u, double v);

/**
 * @brief 计算 B 样条曲线的导数，允许求导次数高于 p 次，这对于 NURBS 曲线是必须的 : ders[j] = C^(j)(u)
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] knots 节点向量
 * @param[in] n 节点向量长度
 * @param[in] p 次数
 * @param[in] d 导数阶数
 * @param[in] u
 * @return EXT_SYMBOL
 */
template <class T>
EXT_SYMBOL std::vector<T> curve_deriv_alg1(const T *cpts, const double *knots, int n, int p, int d, double u);

/**
 * @brief 计算有理 B 样条曲线的导数，允许求导次数高于 p 次，这对于 NURBS 曲线是必须的 : ders[j] = C^(j)(u)
 *
 * @tparam N
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] knots 节点向量
 * @param[in] n 节点向量长度
 * @param[in] p 次数
 * @param[in] d 导数阶数
 * @param[in] u
 * @return EXT_SYMBOL
 */
template <int N, class T>
EXT_SYMBOL std::vector<VecT<N - 1, T>> rat_curve_deriv(const VecT<N, T> *cpts, const double *knots, int n, int p, int d,
                                                       double u);

/**
 * @brief 计算 B 样条曲线的导数，要求 d <= p
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] knots 节点向量
 * @param[in] n 节点向量长度
 * @param[in] p 次数
 * @param[in] d 导数阶数
 * @param[in] u
 * @return EXT_SYMBOL
 */
template <class T>
EXT_SYMBOL std::vector<T> curve_deriv_alg2(const T *cpts, const double *knots, int n, int p, int d, double u);

/**
 * @brief 计算 B 样条曲线导数的控制点，要求 d <= p
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] knots 节点向量
 * @param[in] n 节点向量长度
 * @param[in] p 次数
 * @param[in] d 导数阶数
 * @param[in] r1 起始控制点编号 r1 = 0, r2 = n 时计算所有控制点
 * @param[in] r2 结束控制点编号
 * @return EXT_SYMBOL
 */
template <class T>
EXT_SYMBOL std::vector<std::vector<T>> curve_deriv_cpts(const T *cpts, const double *knots, int n, int p, int d, int r1,
                                                        int r2);

/**
 * @brief 计算 B 样条曲面上的导数，允许求导次数高于 p,q 次，这对于 NURBS 曲面是必须的 : ders[j][k] = S^(j,k)(u,v)
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] uknots u 方向节点向量
 * @param[in] vknots v 方向节点向量
 * @param[in] n u 方向节点向量长度
 * @param[in] m v 方向节点向量长度
 * @param[in] p u 方向次数
 * @param[in] q v 方向次数
 * @param[in] d 导数阶数
 * @param[in] u
 * @param[in] v
 * @return EXT_SYMBOL
 */
template <class T>
EXT_SYMBOL std::vector<std::vector<T>> surface_deriv_alg1(const T *cpts, const double *uknots, const double *vknots,
                                                          int n, int m, int p, int q, int d, double u, double v);

/**
 * @brief 计算有理 B 样条曲面上的导数，允许求导次数高于 p,q 次，这对于 NURBS 曲面是必须的 : ders[j][k] = S^(j,k)(u,v)
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] uknots u 方向节点向量
 * @param[in] vknots v 方向节点向量
 * @param[in] n u 方向节点向量长度
 * @param[in] m v 方向节点向量长度
 * @param[in] p u 方向次数
 * @param[in] q v 方向次数
 * @param[in] d 导数阶数
 * @param[in] u
 * @param[in] v
 * @return EXT_SYMBOL
 */
template <int N, class T>
EXT_SYMBOL std::vector<std::vector<VecT<N - 1, T>>> rat_surface_deriv(const VecT<N, T> *cpts, const double *uknots,
                                                                      const double *vknots, int n, int m, int p, int q,
                                                                      int d, double u, double v);

/**
 * @brief 节点插入算法，可以在任意节点区间中插入节点
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] knots 节点向量
 * @param[in] n 节点向量长度
 * @param[in] p 次数
 * @param[in] u 插入点
 * @param[in] k 插入位置
 * @param[in] s 节点重数
 * @param[in] r 插入次数
 * @param[out] Q 新的控制点
 * @param[out] U 新的节点向量
 */
template <class T>
EXT_SYMBOL void curve_knot_ins(const T *cpts, const double *knots, int n, int p, double u, int k, int s, int r, T *Q,
                               double *U);

/**
 * @brief 节点插入算法，可以在任意节点区间中插入节点
 *
 * @tparam T
 * @tparam D
 * @param[in] cpts 控制点
 * @param[in] knots 节点向量
 * @param[in] n u 方向节点向量长度
 * @param[in] m v 方向节点向量长度
 * @param[in] p u 方向次数
 * @param[in] q v 方向次数
 * @param[in] uv 插入点
 * @param[in] k 插入位置
 * @param[in] s 节点重数
 * @param[in] r 插入次数
 * @param[out] Q 新的控制点
 * @param[out] UV 新的节点向量
 * @return EXT_SYMBOL
 */
template <class T, int D>
EXT_SYMBOL void surface_knot_ins(const T *cpts, const double *knots, int n, int m, int p, int q, double uv, int k,
                                 int s, int r, T *Q, double *UV);

/**
 * @brief 通过节点插入计算 B 样条曲线上的点
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] knots 节点向量
 * @param[in] n 节点向量长度
 * @param[in] p 次数
 * @param[in] u
 * @return T
 */
template <class T> EXT_SYMBOL T curve_point_by_corner_cut(const T *cpts, const double *knots, int n, int p, double u);

/**
 * @brief 曲线节点细化算法，一次插入多个节点
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] knots 节点向量
 * @param[in] n 节点向量长度
 * @param[in] p 次数
 * @param[in] x 插入节点向量
 * @param[in] r 插入节点向量长度
 * @param[out] Q 新的控制点
 * @param[out] U 新的节点向量
 * @return EXT_SYMBOL
 */
template <class T>
EXT_SYMBOL void curve_knot_refine(const T *cpts, const double *knots, int n, int p, const double *x, int r, T *Q,
                                  double *U);

/**
 * @brief 曲面节点细化算法，一次插入多个节点
 *
 * @tparam T
 * @tparam D
 * @param[in] cpts 控制点
 * @param[in] knots 节点向量
 * @param[in] n u 方向节点向量长度
 * @param[in] m v 方向节点向量长度
 * @param[in] p u 方向次数
 * @param[in] q v 方向次数
 * @param[in] x 插入节点向量
 * @param[in] r 插入节点向量长度
 * @param[out] Q 新的控制点
 * @param[out] UV 新的节点向量
 * @return EXT_SYMBOL
 */
template <class T, int D>
EXT_SYMBOL void surface_knot_refine(const T *cpts, const double *knots, int n, int m, int p, int q, const double *x,
                                    int r, T *Q, double *UV);

/**
 * @brief 曲线分解为 Bezier 曲线 : Q[j][k] 表示第 j 段 Bezier 曲线上的第 k 个控制点
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] knots 节点向量
 * @param[in] n 节点向量长度
 * @param[in] p 次数
 * @return EXT_SYMBOL
 */
template <class T>
EXT_SYMBOL std::vector<std::vector<T>> curve_decompose(const T *cpts, const double *knots, int n, int p);

/**
 * @brief 曲线节点删除算法，只有当移除后的误差小于容差时才删除节点
 *
 * @tparam T
 * @tparam Rat
 * @param[in] cpts 控制点
 * @param[in] knots 节点向量
 * @param[in] n 节点向量长度
 * @param[in] p 次数
 * @param[in] u 删除节点
 * @param[in] k 删除位置
 * @param[in] s 节点重数
 * @param[in] r 删除次数
 * @param[int] tol 容差
 * @return EXT_SYMBOL
 */
template <class T, bool Rat>
EXT_SYMBOL int curve_knot_remove(T *cpts, double *knots, int n, int p, double u, int k, int s, int r, double tol = 0.0);

/**
 * @brief 曲线升阶算法，只对于 clamped 类型有效
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] knots 节点向量
 * @param[in] n 节点向量长度
 * @param[in] p 次数
 * @param[in] t 升阶次数
 * @param[out] Q 新的控制点
 * @param[out] U 新的节点向量
 * @param[in] step1 原始控制点步长
 * @param[in] step2 新控制点步长
 * @return EXT_SYMBOL
 */
template <class T>
EXT_SYMBOL int curve_elevate(const T *cpts, const double *knots, int n, int p, int t, T *Q, double *U, int step1 = 1,
                             int step2 = 1);

/**
 * @brief 曲线降阶算法，只对于 clamped 类型有效
 *
 * @tparam T
 * @param[in] cpts 控制点
 * @param[in] knots 节点向量
 * @param[in] n 节点向量长度
 * @param[in] p 次数
 * @param[out] Q 新的控制点
 * @param[out] U 新的节点向量
 * @param[in] tol 容差
 * @return EXT_SYMBOL
 */
template <class T>
EXT_SYMBOL int curve_reduce(const T *cpts, const double *knots, int n, int p, T *Q, double *U, double tol);

} // namespace geom

} // namespace xi
